Phillips  Laboratory  Technical  Report  PL-TR-95-2132  (21  Sep  95) 
Approved  for  Public  Release 

Long-range  Lattice-Gas  Algorithm 


J.  Yepez 

Phillips  Laboratory,  Hanscom  Field,  Massachusetts  01731 

(10  Oct  1994) 


Abstract 

Presented  is  a  novel  algorithmic  method  for  simulating  complex  fluids, 
for  instance  multiphase  single  component  fluids  and  molecular  systems. 
The  algorithm  falls  under  a  class  of  single-instruction  multiple-data  com¬ 
putation  known  as  lattice-gases,  and  therefore  inherits  exact  computabil¬ 
ity  on  a  discrete  spacetime  lattice.  Our  contribution  is  the  use  of  nonlocal 
interactions  that  allow  us  to  model  a  richer  set  of  physical  dynamics,  such 
as  crystallization  processes,  yet  to  do  so  in  a  way  that  remains  locally 
computed.  A  simple  computational  scheme  is  employed  that  allows  all 
the  dynamics  to  be  computed  in  parallel  with  two  additional  bits  of  local 
site  data,  for  outgoing  and  incoming  messengers,  regardless  of  the  number 
of  long-range  neighbors.  The  computational  scheme  is  an  efficient  decom¬ 
position  of  a  lattice-gas  with  many  neighbors.  It  is  conceptually  similar 
to  the  idea  of  virtual  intermediate  particle  momentum  exchanges  that  is 
well  known  in  particle  physics.  All  2-body  interactions  along  a  particular 
direction  define  a  spatial  partition  that  is  updated  in  parallel.  Random 
permutation  through  the  partitions  is  sufficient  to  recover  the  necessary 
isotropy  as  long  as  enough  momentum  exchange  directions  are  used.  The 
algorithm  is  implemented  on  the  CAM-8  prototype. 
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1  Introduction 


To  illustrate  the  how  one  may  implement  a  lattice-gas  with  long-range  interac¬ 
tions,  let  us  consider  for  simplicity  a  two-dimensional  example  with  a  system 
having  only  one  interaction  range  and  consider  only  the  case  of  an  attractive 
interaction.  The  more  general  case  of  multiple  interaction  ranges  with  both  re¬ 
pulsive  and  attractive  interactions  and  in  higher  dimensions  will  follow  directly. 
Specifically,  we  discuss  our  new  molecular  dynamics  lattice-gas  algorithm  that 
uses  eight  interaction  ranges  and  both  repulsive  and  attractive  interactions  to 
approximate  a  Lenord-Jones  intermolecular  potential. 

Our  long-range  lattice-gas  has  been  implemented  on  the  MIT  cellular  au¬ 
tomata  machine  prototype,  the  CAM-8.  Consequently,  we  first  give  a  brief 
description  of  the  CAM-8  architecture.  We  next  give  a  brief  description  of  what 
a  lattice-gas  automaton  is  and  explain  why  it  is  an  exactly  computable  represen¬ 
tation  of  a  dynamical  system.  One  of  the  principle  requirements  for  a  lattice-gas 
with  microscopic  finite-point  group  symmetry  to  give  rise  to  macroscopic  con¬ 
tinuous  rotational  symmetry  is  that  the  underlying  lattice  must  be  isotropic. 
Therefore  we  describe  what  it  means  for  a  lattice  to  be  isotropic.  Working 
in  two-dimensions  is  much  easier  that  working  in  three,  both  for  implement¬ 
ing  computer  models  and  for  describing  them.  For  this  reason  we  present  our 
long-range  lattice-gas  algorithm  in  two-dimensions  on  the  triangular  lattice. 

When  introduced  to  the  triangular  lattice-gas  model  for  the  first  time,  one 
inevitably  asks  the  following  question:  Why  does  the  discrete  dynamics  fail  to 
reproduce  the  correct  continuum  hydrodynamic  limit  when  implemented  on  a 


1 


square  lattice?  One  finds  that  four  momentum  states  are  insufficient  by  noting 
that  the  derivation  of  the  Navier-Stokes  equation  relies  on  the  expansion  of  the 
momentum  flux  density  tensor  in  terms  of  the  isotropic  tensor  E 4  which  itself 
could  be  expanded  in  products  of  two-dimensional  Kronecker  deltas,  given  below 
in  (10).  For  the  square  lattice  case,  the  lattice  vectors  are  orthogonal  and  E4 
cannot  be  decomposed  into  two-dimensional  Kronecker  deltas.  Instead 

Eijki\s=i  =  2(5  ijki 

where  Sijki  is  a  four-dimensional  Kronecker  delta  illustrating  the  lack  of  isotropy 
of  the  momentum  flux  density  on  a  square  lattice-gas.  Since  five  nearest  neigh¬ 
bors  is  not  space  filling,  the  next  possible  choice  is  six  or  the  triangular  lat¬ 
tice.  The  simplest  discrete  dynamics  in  two-dimensions  is  known  as  a  hexag¬ 
onal  lattice-gas  or  an  FHP  lattice-gas,  after  its  originators  Uriel  Frisch,  Brosl 
Hasslacher,  and  Yves  Pomeau  [1].  Since  the  long-range  lattice-gas  still  retains 
local  collisions,  we  present  the  simple  FHP-model  here  for  completeness.  Next 
we  introduce  the  long-range  2-body  interaction,  restricting  ourselves  to  central- 
body  attractive  interactions,  for  the  sake  of  simplicity.  We  discuss  two  different 
bound  states,  the  bounce-back  orbit  and  the  clockwise  orbit. 

When  implementing  lattice-gas  algorithms  it  is  often  useful  try  to  find  eco¬ 
nomical  ways  of  expressing  the  collisions  or  interactions  so  as  to  reduce  the  size  of 
a  look-up  table  or  reduce  the  depth  of  the  logical  representation  of  the  algorithm. 
To  this  end  we  briefly  discuss  some  symmetries  inherent  in  long-range  interac¬ 
tions,  in  particular  we  introduce  partity  and  conjugation  symmetries.  Finally, 
we  discuss  our  implementation  of  a  multi-long-range  lattice-gas.  Remarkably, 
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this  methodology  allows  us  to  model  solid-state  dynamics  and  as  such  offers  an 
alternative  to  the  usual  method  of  computational  molecular  dynamics. 


2  The  Cellular  Automata  Machine  CAM-8 


Figure  1:  MIT  Laboratory  for  Computer  Science  cellular  automata  machine  CAM-8.  This  8 
module  prototype  can  evolve  a  D-dimensional  cellular  space  with  32  million  sites  where  each 
site  has  16  bits  of  data  with  a  site  update  rate  of  200  million  per  second. 


The  cellular  automata  machine  CAM-8  architecture  devised  by  Norman  Mar- 
golus  of  the  MIT  Laboratory  for  Computer  Science  [2,  3]  is  the  latest  in  a  line  of 
cellular  automata  machines  developed  by  the  Information  Mechanics  Group  at 
MIT  [4,  5,  6].  It  is  optimized  for  performing  lattice-gas  simulations.  The  CAM- 
8  architecture  itself  is  a  simple  abstraction  of  lattice  gas  dynamics.  Lattice  gas 
data  streaming  and  collisions  are  directly  implemented  in  the  architecture.  The 


3 


communication  network  is  a  cartesian  three-dimensional  mesh.  Crystallographic 
lattice  geometries  can  be  directly  embedded  into  the  CAM-8.  Each  site  of  the 
lattice  has  a  certain  number  of  bits  (a  multiple  of  16)  which  we  refer  to  as  a 
“cell”.  Each  bit  of  the  cell,  or  equivalently  each  bit  plane  of  the  lattice,  can 
be  translated  through  the  lattice  in  any  arbitrary  direction.  The  translation 
vectors  for  the  bit  planes  are  termed  “kicks”.  The  specification  of  the  x,y,  and 
z  components  of  the  kicks  for  each  bit  plane  (or  hyperplane)  exactly  defines 
the  lattice.  The  kicks  can  be  changed  during  the  simulation.  Thus,  the  data 
movement  in  the  CAM-8  is  general.  Once  the  kicks  are  specified,  the  coding 
of  the  lattice-gas  streaming  is  completed.  In  effect,  the  kicks  determine  all  the 
global  permutations  of  the  data. 

Local  permutations  of  data  occur  within  the  cells.  These  permutations  are 
the  computational  metaphor  for  physical  collisions  between  particles.1  All  local 
permutations  are  implemented  in  look-up  tables.  That  is,  all  possible  physical 
events  with  a  certain  input  configuration  and  a  certain  output  configuration 
are  precomputed  and  stored  in  SRAM,  for  fast  table  look-up.  The  width  of  the 
CAM-8  look-up  tables  are  limited  to  16-bits,  or  64K  entries.  This  is  a  reasonable 
width  satisfying  the  opposing  considerations  of  model  complexity  versus  memory 
size  limitations  for  the  SRAM.  Site  permutations  of  data  wider  than  16-bits 
must  be  implemented  in  several  successive  table  look-up  passes.  Since  the  look¬ 
up  tables  are  double  buffered,  a  scan  of  the  space  can  be  performed  while  a  new 
look-up  table  is  loaded  for  the  next  scan. 

1  Locally,  the  CAM-8  is  not  limited  to  performing  only  permutations,  it  can  do  general 
mappings.  However,  since  we  are  interested  in  only  particle  conserving  reversible  dynamics, 
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Figure  2:  CAM-8  system  diagram,  (a)  A  single  processing  node,  with  DRAM  site  data 
flowing  through  an  SRAM  lookup  table  and  back  into  DRAM,  (b)  Spatial  array  of  CAM-8 
nodes,  with  nearest-neighbor  (mesh)  interconnect  (1  wire/bit-slice  in  each  direction. 


Figure  2  is  a  schematic  diagram  of  a  CAM-8  system.  On  the  left  is  a  sin¬ 
gle  hardware  module — the  elementary  “chunk”  of  the  architecture.  On  the 
right  is  an  indefinitely  extendable  array  of  modules  (drawn  for  convenience  as 
two-dimensional,  the  array  in  normally  three-dimensional).  A  uniform  spatial 
calculation  is  divided  up  evenly  among  these  modules,  with  each  module  sim¬ 
ulating  a  volume  of  up  to  millions  of  fine-grained  spatial  sites  in  a  sequential 
fashion.  In  the  diagram,  the  solid  lines  between  modules  indicate  a  local  mesh 
interconnection.  These  wires  are  used  for  spatial  data  movements.  There  is 
also  a  tree  network  (not  shown)  connecting  all  modules  to  the  front-end  host, 
a  SPARC  workstation  with  a  custom  SBus  interface  card,  controls  the  CAM-8. 
It  downloads  a  bit-mapped  pattern  as  the  initial  condition  for  the  simulations. 

It  also  sends  a  “step-list”  to  the  CAM-8  to  specify  the  sequence  of  kicks  and 
permutations  are  sufficient. 
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scans  that  evolve  the  lattice-gas  in  time.  One  can  view  the  lattice-gas  simu¬ 
lation  in  real-time  since  a  custom  video  module  captures  site  data  for  display 
on  a  VGA  monitor,  a  useful  feature  for  lattice-gas  algorithm  development,  test 
and  evaluation.  The  CAM-8  has  built-in  25-bit  event  counters  allowing  real¬ 
time  measurements  without  slowing  the  lattice-gas  evolution.  This  feature  is 
used  to  do  real-time  coarse-grain  block  averaging  of  the  lattice-gas  number  vari¬ 
ables  and  to  compute  the  components  of  the  momentum  vectors  for  each  block. 

The  amount  of  coarse-grained  data  is  sufficiently  small  to  be  transferred  back 
to  the  front-end  host  for  graphical  display  as  an  evolving  flow  field  within  an 
X- window. 

3  Lattice-Gas  Automaton:  An  Exactly  Computable 
Dynamical  System 

A  boolean  formulation  of  an  exactly  computable  dynamical  system,  known  as 
a  lattice-gas,  may  be  stated  in  a  way  that  is  consistent  with  the  Boltzmann 
equation  for  kinetic  transport.  In  essence  the  lattice-gas  dynamics  are  a  sim¬ 
plified  form  of  molecular  transport  as  we  restrict  ourselves  to  a  cellular  phase 
space.  The  macroscopic  equations,  in  particular  the  continuity  equation  and 
the  Navier-Stokes  equation,  are  obtained  by  coarse-graining  over  a  discrete  mi- 
crodynamical  transport  equation  for  number  boolean  variables.  The  scheme 
employs  the  finite-point  group  symmetry  of  a  crystallographic  spatial  lattice.  It 
is  somewhat  inevitable  that  to  obtain  an  exactly  computable  representation  of 
fluid  dynamics  one  must  perform  a  statistical  treatment  over  discrete  number 
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variables. 


Before  introducing  the  basic  lattice-gas  microdynamical  transport  equation, 
let  us  give  some  notational  conventions.  We  consider  a  spatial  lattice  with  N 
total  sites.  The  fundamental  unit  of  length  is  the  size  of  a  lattice  cell,  l,  and 
the  fundamental  unit  of  time,  r,  is  the  time  it  takes  for  a  speed-one  particle  to 
go  from  one  lattice  site  to  a  nearest  neighboring  site.  Particles,  with  unit  mass 
to,  propagate  on  the  lattice.  The  unit  lattice  propagation  speed  is  denoted 
by  c  =  Particles  occupy  this  discrete  space  and  can  have  only  a  finite  B 
number  of  possible  momenta.  The  lattice  vectors  are  denoted  by  eai  where 
a  =  1,2 For  example,  for  a  single-speed  gas  on  a  triangular  lattice, 
a  =  1,2,...,  6.  A  particle’s  state  is  completely  specified  at  some  time,  t,  by 
specifying  its  position  on  the  lattice,  Xi,  and  its  momentum,  pi  =  mceai,  at  that 
position.  The  particles  obey  Pauli  exclusion  since  only  one  particle  can  occupy 
a  single  momentum  state  at  a  time.  The  total  number  of  configurations  per  site 
is  2s.  The  total  number  of  possible  single  particle  momentum  states  available 
in  the  system  is  lVtotal  =  BN.  With  P  particles  in  the  system,  we  denote  the 
filling  fraction  by  d  =  ^ . 

The  number  variable,  denoted  by  na(x,  t),  takes  the  value  of  one  if  a  particle 
exists  at  site  x  at  time  t  in  momentum  state  m.cea,  and  takes  the  value  of  zero 
otherwise.  The  evolution  of  the  lattice-gas  can  then  be  written  in  terms  of  na  as 
a  two-part  process:  a  collision  and  streaming  part.  The  collision  part  reorders 
the  particles  locally  at  each  site. 

«'a(x,i)  =  na(x,t)  +  f2a(n(x,t)),  (1) 
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where  f la  represents  the  collision  operator  and  in  general  depends  on  all  the 
particles,  n  at  the  site.  So  as  a  short-hand  we  suppress  the  index  on  the  occu¬ 
pation  variable  when  it  is  an  argument  of  fla(n(x,t))  to  represent  this  general 
dependence.  In  the  streaming  part  of  the  evolution  the  particle  at  position  x 
“hops”  to  its  neighboring  site  at  x  +  lea  and  then  time  is  incremented  by  r 

n'a(x  +  lea,  t  +  t)  =  nQ(x,  t)  +  (2a(n(x,  t)).  (2) 

Equation  (2)  is  the  lattice-gas  microdynamical  transport  equation  of  motion. 

The  collision  operator  can  only  permute  the  particles  locally  on  the  site  since 
we  wish  the  local  particle  number  to  be  conserved  before  and  after  the  collision. 
That  is, 

n(x,t)  =  ^Va(x,t)  =  ^na(x,f).  (3) 

a  a 

Equation  (3)  defines  the  local  number  density.  Summing  (2)  over  lattice  direc¬ 
tions  then  implies  the  following  constraint  on  the  collision  operator 

5>»  =  °-  (4) 

a 

We  may  define  the  local  momentum  as 

Pi (x,  t)  =  mc^eain'0(x,f)  =  mc^  ea»n0(x,  t),  (5) 

a  a 

which  of  course  must  also  be  conserved  before  and  after  a  collision.  Again,  this 
imposes  a  constraint  on  the  collision  operator. 

]Teafia  =  0.  (6) 

a 

As  a  matter  of  notation  it  should  be  understood  that  whenever  a  directionally 
dependent  quantity  is  written,  its  subscripted  index  is  taken  modulo  B.  Using 


the  number  variable,  for  example,  it  is  understood  that 


^a+6  ^  mod  b  (a+6)  • 


(7) 


As  a  short  hand,  a  negative  indice  will  represent  the  antiparallel  direction,  so 
since  ea+B  =  — ea  we  may  write 


4  Isotropic  Lattice  Tensors 


(8) 


We  construct  an  n-th  rank  tensor  composed  of  a  product  of  lattice  vectors  [7] 
£(n)  =  Eix...in  =  •  •  •  (e„)iB,  (9) 

a 

where  a  =  1, . . . ,  B.  All  odd  rank  E  vanish.  We  wish  to  express  E^2n'>  in  terms 
of  Kronecker  deltas,  6tj  =  1  for  i  =  j  and  zero  otherwise.  We  can  turn  this 
problem  of  expressing  the  .E-tensors  in  terms  of  products  of  Kronecker  deltas 
into  a  problem  of  combinatoric  counting.  We  use  the  following  tensors 


A2- 

tj 

II 

Oo 

(10) 

A4 

^ ijkl 

—  $ij$kl  “1“  ^ik^jl  “1“  ^il^kj 

(ii) 

and  so  forth.  Then  we  know  that  if  E  is  isotropic  it  must  be  proportional  to  A 

E(2n)  oc  A(2n)  (12) 


and  that  the  constant  of  proportionality  may  be  obtained  by  counting  the  num¬ 
ber  of  ways  we  could  write  a  term  comprising  a  product  of  n  Knonecker  deltas. 
Consider  for  example  the  case  n  =  2.  Since  the  Knonecker  delta  is  symmetric  in 
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its  indices,  the  following  four  products  are  identical:  SijSki  =  SijSik  =  SjjSki  = 
fijidik-  The  degeneracy  is  22.  Furthermore,  the  order  of  the  Kronecker  deltas 
also  doesn’t  matter  since  they  commute;  that  is,  SijSki  =  The  degen¬ 

eracy  is  2!.  For  the  case  where  n  is  arbitrary,  there  are  2n  identical  ways  of 
writing  the  product  of  n  Kronecker  deltas.  For  each  choice  of  indices,  there  are 
an  additional  n!  number  of  ways  of  ordering  the  products.  Therefore,  the  total 
number  of  degeneracies  equals  2 nn!  =  (2 n)!!.  The  total  number  of  permutations 
for  2 n  indices  equals  (2n)!.  So  from  this  counting  procedure  we  know  that  A(2n') 
consists  of  a  sum  of  ^2y.<  =  —  1)!!  terms. 

The  following  relations  will  be  very  useful  throughout  later  developments 


El 

=  0 

(13) 

E2 

B, 

~ 

(14) 

E3 

=  0 

(15) 

E4 

—  _|_  2)  T  Mu  T 

(16) 

In  general,  the  lattice  tensors  are 


E2n+1  =  0 

^  2n  _ _ _ A2n 

D(D  +  2)  •  •  •  (D  +  2n  -  2) 

5  Triangular  Lattice 


(17) 

(18) 


In  a  triangular  lattice  there  are  six  vectors  that  we  enumerate  in  this  section  by 
the  following  convention 


„  /  .  7ra  7 ro 

ea  =  —  (J3in  —  ,cos  — 
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where  a  =  1,2,. ..,6.  The  spatial  coordinates  of  the  lattice  sites  may  be  ex- 


(4,0) 

(4,2) 

1 

; 

•--(43)- 

(3,0) 

| 

(3,2) 

a  =  6 

; 

.(3,1).- 

'-.(3,3). 

a  =  5 

CJL 

a  =  1 

(2,0) 

| 

(2,2) 

! 

a  =  4 

'T- 

ft! 

II 

(1,0) 

T' 

(1,2) 

M2,3) 

r 

a  =  3 

--XU). 

(0,0) 

(0,2) 

i 

-ip,]).- 

j 

-•J0;3i 

(4,4) 


(3,4) 


(2,4) 


(1,4) 


(0,4) 


V3/_ 
2  1 


(a) 


(b) 


Figure  3:  (a)  Lattice  vector  label  convention;  (b)  Hexagonal  lattice  convention  with  lattice 
directions  a  =  3  up  and  a  =  6  down.  Coordinates  above  the  lattice  nodes  are  ( i,j )  memory 
array  indices. 


pressed  as  follows 


(Vs  i 

x*i  =  oU  mod  2) 


(20) 


2  •”  2 

where  i  and  j  are  rectilinear  indices  which  specify  the  data  memory  array  loca¬ 
tion  used  to  store  the  lattice-gas  site  data. 

Let  s  =  ( j  mod  2)(r  mod  2).  Given  a  particle  at  site  it  may  be  shifted 
to  a  site  r  lattice  units  away  to  a  remote  site  ( i',j' )  by  the  following  mapping 


(<',/)  1  = 

(*+  2  -) 

(21) 

(i'J'h  = 

(22) 

3  = 

(*  -  r,j) 

(23) 

= 

( .  r  + 1  .  \ 

l*+  2  s,j+rl 

(24) 
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= 

(»-^-»J  +  r) 

(25) 

{i  +  r,j) 

(26) 

where  ( denotes  the  shifted  site,  i.e.  ( i,j )  —>  with  a  shift  along 

vector  r  =  rea  and  where  division  by  2  is  considered  integer  division. 

6  Local  Collision  Rules 


In  two-dimensions  we  may  use  a  triangular  lattice,  with  six  bits  per  site  encoding 
the  occupation  numbers  of  the  six  possible  momentum  states.  Let  na  be  the 
input  bits  and  n'a  be  the  output  bits  of  a  local  collision.  A  general  collision 
operator  is  constructed  as  follows 


fia  =  5>Qa((CJ),  (27) 

{Ci} 

where  {£;}  is  a  set  of  occupied  particle  states  and  a  =  ±1  is  a  scalar  coefficient 


and  where  each  term  in  the  sum  is  written  in  factorized  form  as 

B 


Qa(il,...,ik)  =  I^a+il 


na+ik 


n(!-n8+j).  (28) 

^o+ii  1  ^a+ik  -  1 


Then  the  FHP  collision  operator  is  the  following: 


«aHP  =  \Qa{  1, 4)  +  ^Qa( 2, 5)  -  Qa( 0, 3)  +  Qa(  1, 3, 5)  -  Qa{ o,  2, 4)  (29) 

n0  =  ^mn4(l  -  n0)(l  -  n2)(l  -  n3)(l  -  n5)  + 

^n2n5(  1  -  n0)(l  -  ni)(l  -  n3)(l  -  n4)  - 
n0n3(l  -  ni)(l  -  n2)(l  -  n4)(l  -  n5)  + 


nin3n5(l  -  n0)(  1  -  n2)(l  -  n4)  - 


non2n4(l  -  ni)(l  -  n3)(l  -  n5) 
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Table  1:  Simple  right-handed  collision  table 


n0 

ni 

n2 

ri3 

77-4 

n5 

no 

n[ 

n'2 

n3 

n'A 

< 

1 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

1 

0 

1 

0 

0 

l 

0 

1 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

1 

0 

1 

0 

0 

1 

0 

1 

0 

1 

0 

i 

0 

0 

1 

0 

1 

0 

1 

0 

1 

0 

1 

0 

1 

1 

0 

1 

0 

1 

0 

Table  2:  Simple  left-handed  collision  table 


n0 

Til 

n2 

n3 

77-4 

n5 

no 

ri'i 

n'2 

n3 

n'4 

«5 

1 

0 

0 

1 

0 

0 

0 

1 

0 

0 

1 

0 

0 

1 

0 

0 

l 

0 

0 

0 

1 

0 

0 

1 

0 

0 

1 

0 

0 

1 

1 

0 

0 

1 

0 

0 

1 

0 

1 

0 

i 

0 

0 

1 

0 

1 

0 

1 

0 

1 

0 

1 

0 

1 

1 

0 

1 

0 

1 

0 

Note  that  it  is  sufficient  to  give  only  fio  since  the  other  components  of  the 
collision  operator  can  be  obtained  simply  by  incrementing  the  indices  of  Qq 
owing  to  the  six-fold  symmetry  of  the  collisions.  The  factors  of  one-half  in  (29) 
are  transition  probabilities  for  the  2-body  collisions,  indicating  a  coin  toss  is 
performed  to  choose  between  even  or  odd  chirality. 

The  possible  two-body  and  three-body  collisions  represented  by  (29)  are 
illustrated  in  figure  (4).  For  two-dimensional  flow,  there  are  four  invariants, 
the  mass,  two  components  of  the  momentum,  and  the  energy.  With  only  the  2- 
body  collision  in  figure  (4),  there  is  an  additional  invariant:  the  difference  in  the 
particle  number  along  each  of  the  three  lattice  directions.  The  3-body  collisions 
in  figure  (4)  were  include  in  the  FHP-model  to  remove  this  spurious  invariant. 
Consequently,  the  collisions  enumerated  in  figure  (4)  are  the  minimally  sufficient 
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100100 ->001001  t 

A  f 

010010 ->  100100  ^ 

{  0  0  1  00  1  ->010010  \ 

Even  Chirality 

100100 ->  010010 

^  010010 ->  001001  } 

,  t 

;  001001  ->  100100 

t  "" 

Odd  Chirality 

1  101010 ->010101  'V'' 

010101 ->101010  J 

Figure  4:  Enumeration  of  FHP  2-body  collisions,  even  and  odd  chirality,  and  3-body  colli¬ 
sions. 


set  to  produce  hydrodynamic  behavior  in  the  continuum  limit. 


7  Long-Range  2-Body  Interactions 


Figure  5:  Simple  bound-state  orbits  due  to  a  long-range  attractive  interaction  where  the 
dotted  path  indicates  the  particle’s  closed  trajectory:  (a)  partition  directions;  (b)  bounce-back 
orbit  with  |Ap|  =  2  and  zero  angular  momentum;  and  (c)  clockwise  orbit  with  |Ap|  =  1  and 
one  unit  of  angular  momentum.  Head  of  the  dashed  arrows  indicates  particles  entering  the 
sites  of  partition  ro  at  time  t.  Tail  of  the  black  arrows  indicates  particles  leaving  those  sites 
at  time  t  +  r.  The  counter-clockwise  orbit  is  not  shown. 


An  interparticle  potential,  F(x  —  x'),  acts  on  particles  spatially  separated 
by  a  fixed  distance,  x  x'  =  r.  An  effective  interparticle  force  is  caused  by  a 
non-local  exchange  of  momentum.  Momentum  conservation  is  violated  locally, 
yet  it  is  exactly  conserved  in  the  global  dynamics. 
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Table  3:  Lattice  vector  components 


a 

x-component 

y-component 

0 

-1 

0 

1 

_  1 

a/3 

2 

2 

1 

2 

VS 

2 

2 

3 

1 

0 

4 

1 

_VS 

K 

2 

1 

2 

V3 

o 

2 

2 

For  the  case  of  an  attractive  interaction,  there  exists  a  bound  states  in  which 
two  particles  orbit  one  another.  Since  the  particle  dynamics  are  constrained  by 
a  crystallographic  lattice  we  expect  polygonal  orbits.  In  figure  5  we  have  de¬ 
picted  two  such  orbits  for  a  hexagonal  lattice-gas.  The  range  of  the  interaction  is 
r.  Two-body  single  range  attractive  interactions  are  depicted  in  figures  5b  and 
5c,  the  bounce-back  and  clockwise  orbits  respectively.  Momentum  exchanges 
occur  along  the  principle  directions.  The  interaction  potential  is  not  spherically 
symmetric,  but  has  an  angular  anisotropy.  In  general,  it  acts  only  on  a  finite 
number  of  points  on  a  shell  of  radius  | .  The  number  of  lattice  partitions  neces¬ 
sary  per  site  is  half  the  lattice  coordination  number,  since  two  particles  lie  on  a 
line.  Though  microscopically  the  potential  is  anisotropic,  in  the  continuum  limit 
obtained  after  coarse-grain  averaging,  numerical  simulation  indicates  isotropy 
is  recovered. 
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8  A  Simple  Example:  Bounce-Back  Orbit 


A  long-range  lattice-gas  of  the  type  we  are  considering  still  possesses  the  usual 
local  dynamics  of  a  hydrodynamic  lattice-gas.  To  extend  the  local  lattice-gas 
update  rules  to  include  long-range  interactions,  we  use  two  additional  bits  of 
local  site  data.  This  will  allow  us  to  implement  a  long-range  interaction  using 
strictly  local  updating  and  therefore  our  algorithm  will  still  be  parallelizable 
just  as  a  usual  local  lattice-gas.  The  two  additional  bits  will  denote  the  occupa¬ 
tion  numbers  of  messenger  particles,  or  “photons” .  The  idea  of  using  messenger 
particles  was  introduced  by  Appert  et  al.  [8] .  We  have  two  types  of  messen¬ 
ger  states,  to  represent  incoming  and  outgoing  conditions,  and  we  denote  the 
messengers  as  zi  and  zr. 

Now  for  the  simplest  long-range  lattice-gas  model,  we  therefore  use  eight  bits 
of  local  site  data.  Since  long-range  interactions  occur  between  remote  spatial 
sites,  say  x  and  x! ,  the  messenger  particles  will  travel  either  parallel  or  antipar¬ 
allel  to  the  vector  r  =  x  —  a?.  All  pairs  of  sites  throughout  the  entire  space  that 
are  separated  by  vector  r  can  therefore  all  be  updated  in  parallel.  We  refer  an 
update  step  of  all  pairs  of  2-body  interactions  along  direction  r  as  a  partition, 
denoted  by  Tr.  All  possible  two-body  interaction  pairs  are  then  computed  by 
performing  all  possible  partitions  of  the  space.  So  it  requires  many  scans  for 
the  space  to  perform  a  single  long-range  interaction  step. 

In  our  two-dimensional  example  using  a  triangular  lattice,  there  are  three 
possible  partitions.  The  number  of  partitions  is  never  smaller  than  half  the  lat- 
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tice  coordination  number.  In  the  two-dimensional  case,  the  simplest  long-range 
lattice-gas  algorithm,  though  perhaps  not  the  most  efficient  algorithm,  is  to  use 
three  sequential  scans  of  the  space,  each  scan  performing  the  updating  necessary 
for  a  single  partition,  see  figure  5a.  Often  times,  depending  on  the  complexity 
of  the  long-range  interactions  and  the  dimensionality  of  the  lattice,  it  is  possible 
to  perform  simultaneous  updating  of  multiple  partitions.  This  of  course  is  more 
desirable  yet  causes  more  complexity.  Furthermore,  this  updating  requires  an 
extra  pair  of  messenger  particles  for  each  partition  to  be  simultaneously  updated. 
For  simplicity,  we  will  not  deal  with  this  case  here,  however  our  implementa¬ 
tion  on  the  CAM-8  does  use  simultaneous  partition  updating — repulsive  and 
attractive  partitions  are  performed  in  parallel  using  four  messenger  bits. 

Let  us  consider  a  simple  example  of  the  long-range  lattice-gas  algorithm, 
the  minimal  model  of  Appert.  Here  we  consider  only  bounce-back  attractive 
interactions.  Suppose  there  is  a  single  particle  at  site  x  =  0  and  there  is  also 
a  single  particle  at  site  x'  =  ri;  that  is,  no(x)  =  1,  n 3(2)  =  0,  n 0(2?)  =  0  and 
713(2?')  =  1  with  all  other  bits  at  x  and  2?'  being  zero,  see  figure  5b.  Here  we  are 
using  the  bit  convention  shown  in  table  3.  Then  the  two  particles  are  separated 
by  a  distance  r  and  are  moving  away  from  each  other.  The  attractive  long-range 
interaction  will  effectively  flip  their  respective  directions  making  710(2?)  =  0, 
713(0;)  =  1,  710(2?')  =  1  and  713(2?')  =  0  so  that  the  two  particles  will  now  be 
moving  toward  each  other.  There  is  a  local  momentum  change  of  2mci  at  x1 
and  an  opposite  momentum  change  of  — 2mci  at  x.  Locally  momentum  is  not 
conserved,  but  nonlocally  it  is. 
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The  first  step  of  the  long-range  interaction  sequence  is  to  choose  a  partition, 
say  Tr,  and  then  to  emit  messenger  particles  along  the  partition  axis.  The  basic 
local  rule  for  this  first  step  is  the  following:  a  photon  is  emitted  at  a  site  if 
there  exists  a  particle  at  that  site  that  can  partake  in  a  long-range  interaction. 
Another  way  of  expressing  this  rule  is:  send  only  if  you  can  receive.  Obviously, 
for  a  particle  to  partake  in  an  interaction  there  must  be  both  a  particle  and  a 
hole  at  that  site.  The  factorized  probability  of  having  such  a  situation  is  just 
d(l  —  d).  So  to  continue  with  our  example,  for  a  photon  to  be  emitted  at  some 
site  x  parallel  or  antiparallel  to  a  partition  direction  i,  we  use  the  following  rule 


Zr(x)  =  n0(x)(l  -  n3(x)) 

(30) 

zi(x)  =  n3(x)(l  -  n0(x)). 

(31) 

Note  that  according  to  this  local  rule,  only  one  photon  can  be  created  at  a  site, 
and  consequently  we  eliminate  the  possibility  of  a  long-range  interaction,  say  of 
range  2 r,  mediated  through  a  doubly  occupied  site.  The  important  consequence 
of  the  emission  step,  is  that  for  two  sites  separated  by  the  interaction  distance, 
r,  if  both  sites  send  photons,  both  will  necessarily  receive  them,  which  strictly 
enforces  nonlocal  momentum  conservation.  Give  and  ye  shall  receive  (provided 
your’s  is  received).  Letting  za  =  zr  and  Z-a  =  Zi,  in  general  we  can  write  the 
emission  step  of  the  minimal  interaction  as 

za(x)  =  n-a(x)(l  -  na(x)),  (32) 

where  a  =  0,1,2  covers  all  the  partitions. 

After  the  emission  step,  follows  a  long-range  kick  of  the  messenger  bits.  In 
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the  simple  example,  all  photons  Zi  are  kicked  along  — ri  and  all  photons  zr  are 
kicked  along  ri.  In  general  for  the  long-range  kick  we  have 

z'a(x  +  rea)  =  za(  x).  (33) 

Finally,  we  have  the  absorption  step  of  the  long-range  interaction  sequence. 
Here  the  local  particle  momentum  state  is  updated  as  the  particles  flip  their 
directions  in  our  example 

n'3(x)  =  n3(x)  +  z((f)n0(f)(l  -  n3(x))  -  z'r(x)n3(x)(l  -  n0(f)) 

(34) 

n'0(x)  =  n0(x)  +  z'r{x)n3(x)(l  -  n0(x))  -  ^(f)?i0(f)(l  -  n3(f)). 

(35) 

Moreover,  in  this  step  all  the  messenger  bits  are  set  to  zero  throughout  the 
entire  space.  For  any  direction,  the  local  absorption  rule  could  more  simply  be 
written  as 

n'a(x)  =  na(x )  +  z'_a(x)za(x )  -  z'a(x)z-a(x).  (36) 

Substituting  in  (32)  and  (33)  into  (36),  we  have  a  single  boolean  expression  in 
terms  of  number  variables  for  a  single  long-range  interaction  step  for  partition 
rr  as  follows 

n'a{x)  =  na(x)  + 

na(x  +  re0)(  1  -  n_a(f  +  rea))n_a(£)(l  -  na(x ))  - 
n-a(x-  rea)(l  -  na(x  -  rea))na(x)(  1  -  n_a(£)) 

(37) 
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Table  4:  Long-range  interaction  sequence 


events 

na(x) 

zi(x) 

zr(x) 

na(x') 

zi(x') 

zr(x') 

initial 

100000 

0 

0 

000100 

0 

0 

emit 

100000 

0 

1 

000100 

1 

0 

kick 

100000 

1 

0 

000100 

0 

1 

absorb 

000100 

0 

0 

100000 

0 

0 

For  convenience  we  define  a  long-range  collision  operator,  Pa,  as  follows 


Pa(x)  =  z'_a(x)za(x), 


(38) 


so  that  we  may  write 

n'a(x)  =  na(x)  +  Pa(x)  -  P-a{x). 


(39) 


The  state  data  for  this  simple  example  we  have  been  considering  is  given  in 
table  4,  which  represents  all  the  steps  of  a  long-range  interaction  sequence  for 
a  partition  along  the  x-axis. 

9  Another  Example:  Clockwise  Orbit 


To  continue  illustrating  our  implementation  of  a  long-range  lattice-gas,  in  this 
section  we  again  consider  a  system  with  a  single  attractive  interaction  of  range 
r,  however  the  local  momentum  states  participating  in  the  interaction  are  not 
along  the  partition  direction.  Yet  in  the  example  given  here,  the  momentum 
exchange  is  still  along  the  partition  direction  so  that  the  long-range  interac¬ 
tion  remains  a  central-body  one,  resulting  in  a  bound  state  with  two  particles 
trapped  in  a  clockwise  orbit.  (Note  that  the  restriction  to  central-body  forces 
is  not  necessary,  but  is  presented  here  for  convenience.)  In  this  slightly  more 
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complicated  example,  the  local  rules  for  photon  emission,  and  absorption,  (32) 
and  (36)  respectively,  have  a  more  general  form  with  the  implication  that  the 
emission  and  absorption  of  photons  is  different  from  the  previous  example  of 
the  bounce-back  orbit  and  should  be  noted  when  making  look-up  tables  to  do 
this  computation.  The  local  photon  emission  rules  can  be  written 


Za(x)  = 

nc{x)(  1  -  nd{x)) 

(40) 

-a{x) 

ng{ £)(1  -  nh(x)) 

(41) 

where  the  bits  c,d,  g,h  must  by  chosen  so  momentum  is  conserved 

ec  -  ed  +  eg  —  eh  =  0  (42) 

as  well  as  be  constrained  by  central-body  parallel  and  perpendicular  momentum 
exchange  conditions 


(ec  -  ed  -  eg  +  eh)  ■  r  =  2A p  (43) 

(ec  -  ed  -  eg  +  eh)  x  f  =  0,  (44) 

where  A p  is  the  momentum  change  per  site  due  to  the  long-range  interaction. 
In  (40)  and  (41)  the  difference,  over  our  previous  example  of  the  bounce-back 
orbit,  is  the  possibility  of  having  two-photons  emitted  at  a  single  site. 

To  be  explicit,  for  the  two-dimensional  triangular  lattice,  we  can  satisfy  (42), 
(43),  and  (44)  by  choosing  the  indices  c,  d ,  g ,  h  as  follows: 


a  —  2 

(45) 

a  —  1 

(46) 

21 


9 


(47) 


=  — c 

h  =  —d.  (48) 

An  example  of  this  choice  of  indices  is  illustrated  in  figure  5c.  Then  the  emission 
rule,  (40)  and  (41),  is  simply 

za(x)  =  na-2(x)(  1  -  na-i(x))  (49) 

Since  the  kicking  of  the  photons  is  the  same  in  this  example  as  in  the  previous 
one,  (33)  still  holds 

z'a(x  +  rea)  =  za(x). 

By  re-expressing  (36)  more  generally,  we  can  write  a  local  absorption  rule 

n'a(x)  =  na(x)  +  Z'_{a+1)(x)za+1(x)  -  <_i(:r)z_(a- i)(z)  (50) 

or  more  elegantly 

n'a(x)  =  na(x)  +  Pa+i(x)  -  P-a+i{x).  (51) 

Substituting  in  (49)  and  (33)  into  (50)  and  after  some  manipulation  of  the 
indices,  we  have  a  single  boolean  expression  in  terms  of  number  variables  for  a 
single  long-range  interaction  step  for  partition  Tr  as  follows 

n'a{x)  =  na(x)  + 

na+2(x  +  rea+i)(l  -  n-a(x  +  rea+i))na_i(f)(l  -  na(x))  - 
n-a(x  -  rea- i)(l  -  na_2(f  -  rea_i))n0(£)(l  -  na+i(£)). 

Table  5  gives  the  local  site  data  for  the  x-axis  partition  of  a  clock- wise  orbit. 
The  particle  n±{x)  acts  as  a  kind  of  spectator  is  this  example,  illustrating  that 
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Table  5:  Long-range  interaction  sequence  with  two  photons  emitted  at  a  single 
site 


labels 

na(x) 

zi(x) 

zr(x) 

na{x') 

zi{x') 

zr(x') 

events 

010010 

0 

0 

000010 

0 

0 

emit 

010010 

1 

1 

000010 

1 

0 

kick 

010010 

1 

0 

000010 

0 

1 

absorb 

001010 

0 

0 

000001 

0 

0 

Table  6:  Long-range  interaction  sequence  with  two  photons  emitted  and  ab¬ 
sorbed  at  site  x'  in  a  back-to-back  interaction 


events 

na(x) 

zi{x) 

zr(x) 

na(x') 

zi(x') 

zr(x') 

na(x”) 

zi(x") 

zr(x ") 

initial 

010000 

0 

0 

010010 

0 

0 

000010 

0 

0 

emit 

010000 

0 

1 

010010 

1 

1 

000010 

1 

0 

kick 

010000 

1 

0 

010010 

1 

1 

000010 

0 

1 

absorb 

001000 

0 

0 

001001 

0 

0 

000001 

0 

0 

two  photons  can  be  emitted  from  a  single  site.  It  is  also  possible  to  have  two 
photons  absorbed  at  a  single  site.  Let  us  consider  a  back-to-back  interaction  over 
three  sites.  Suppose  there  are  particles  at  sites  x  =  0,  x'  =  ri,  and  x"  =  2ri. 
Table  6  gives  the  site  data  for  these  sites  were  there  are  two  photons  emitted 
and  absorbed  at  x'  in  the  middle  location. 

The  minimal  model  using  only  a  attractive  interactions  models  a  fluid  with 
liquid  and  gas  phases  at  zero  temperature.  Fiqure  6  shows  the  time  evolution  of 
the  phase  separation  process  in  this  case  at  a  density  d  =  0.07  and  interaction 
range  r  =  61  and  illustrates  the  type  of  physical  simulation  that  can  be  achieved 
with  the  simplest  long-range  lattice-gas  algorithm. 
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Figure  6:  Time  evolution  of  a  liquid-gas  phase  separation  for  a  lattice-gas  with  long  range 
attractive  interactions  at  range  r  =  6  on  a  1024  x  1024  lattice  starting  with  a  uniformly 
random  configuration  of  density  d  =  0.07. 
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Figure  7:  Examples  of  two-body  finite  impact  parameter  collisions  along  the 
ro-direction.  The  four  terms  of  interaction  Hamiltonian  for  |Ap|  =  1.  Dashed 
arrows  depict  input  states  and  black  arrows  depict  output  states. 
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10  Symmetries  in  Long  Range  Interactions 


For  the  case  of  an  attractive  interaction  as  mentioned  earlier,  we  would  expect 
that  there  exists  a  bound  state  in  which  two  particles  orbit  one  another,  for 
example  the  clockwise  orbit  shown  in  figure  5c.  The  range  of  the  interaction  is 
r.  This  hexagonal  orbit  is  also  shown  in  the  bottom-right  corner  of  figure  7.  In 
this  figure  the  hexagon’s  radius  is  also  labeled  as  r  and  should  not  be  confused 
with  the  interaction  range  which  is  the  hexagon’s  diameter.  Simlarly  a  counter¬ 
clockwise  orbit  is  possible,  see  the  top-left  corner  of  figure  7.  A  time-reversal 
invariance  exists  between  these  two  cases  with  respect  to  conjugation  and  parity 
operations  as  depicted  in  figure  7. 

These  diagrams  describe  the  possible  2-body  collisions  that  can  be  so  gener¬ 
ated,  including  repulsive  interactions.  This  logical  correspondence  between  the 
different  type  of  2-body  interactions,  allows  one  to  achieve  a  more  efficient  im¬ 
plementation  of  a  long-range  lattice-gas  algorithm  than  what  we  have  achieved 
on  the  CAM-8  since  we  have  not  used  this  form  of  logical  economy.  Further¬ 
more,  there  is  a  correspondance  property  for  r  =  0,  the  situation  reduces  to  the 
2-body  collisions  in  the  FHP  lattice-gas  2 . 

The  long  range  interactions  considered  here  have  the  following  properties 
which  simplify  a  computational  implementation:  1)  there  exists  only  parallel 
momentum  exchange,  a  restriction  for  modelling  central  body  forces;  2)  the 

interaction  acts  only  along  the  principal  lattice  directions;  3)  there  exists  time- 

2For  r  =  0,  the  |Ap|  =  1  interactions  reduce  to  a  rotation  of  the  states,  7t.(  ^ )  and  'R,(  -7 ) : 
and  the  |Ap|  =  2  interactions  reduce  to  the  identity  operation. 
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reversal  invariance  with  respect  a  conjugation  and  parity  operations;  and  4)  a 
bias  to  the  interactions  can  be  assigned  by  coupling  to  heat-bath  states.  In  our 


previous  examples,  the  bounce-back  and  clockwise  orbits,  we  discussed  proper¬ 
ties  1  and  2.  In  this  section  we  have  discussed  property  3.  We  have  discussed 
property  4  elsewhere  [9] . 


11  Central-Body  Interaction  Neighborhood  and 
2D  Crystallization 


(-R ,  -2R) 


Figure  8:  Twelve  neighbors  on  a  triangular  lattice  that  can  participate  in  long-range  central- 

body  interactions.  Interaction  range  D  =  7  is  depicted,  where  R  =  ~  4  to  within 

v  3 

1.03  percent  error.  Computer  memory  space  coordinates  ( x,y )  are  given  adjacent  to  each 
neighboring  site. 


In  the  previous  two  sections  we  illustrated  two  examples  of  the  long-range  in- 
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teractions,  the  bounce-back  and  clockwise  orbit  on  a  two-dimensional  triangular 
lattice.  In  general,  the  long-range  interaction  step  will  involve  many  partitions, 
both  attractive  and  repulsive  interactions,  and  multiple  ranges.  In  our  CAM-8 
implementation  of  a  long-range  lattice-gas  with  central-body  interactions,  we 
actually  use  12  neighbors  in  two-dimensions,  see  figure  8.  The  triangular  lat¬ 
tice  is  superposed  over  a  square  lattice,  which  appears  rhomboidal  in  the  figure. 
The  square  lattice  is  often  used  for  embedding  the  site  data  into  computer  mem¬ 
ory,  which  is  rectilinear.  This  kind  of  embedding  is  the  simplest  and  used  for 
simulations  that  possess  periodic  boundary  conditions.  The  reason  for  using  12 
neighbors  is  to  try  to  achieve  a  higher  degree  of  local  symmetry.  In  doing  molec¬ 
ular  dynamics  modeling  with  a  multi-long-range  lattice-gas,  we  have  found  that 
12  neighbors  are  necessary  to  recover  macroscopic  isotropy.  In  particular,  12 
neighbors  are  necessary  to  have  the  emergent  crystalline  solid  be  able  to  freely 
rotate  in  space.  A  mean-field  analysis  of  the  latttice-gas  crystallization  method 
has  been  presented  elsewhere  [10].  Figure  8  shows  a  ring  of  range  7  lattice 
spacings. 

To  implement  the  crystallization  algorithm,  we  use  up  to  eight  ranges  in 
two-dimensions.  That  is,  eight  rings  of  the  type  shown  in  figure  8  for  a  total  of 
96  neighbors.  Half  the  rings  are  used  for  attractive  interactions  and  the  other 
half  are  used  for  repulsive  interactions.  Typically,  the  inner  rings  are  attrac¬ 
tive,  the  middle  rings  are  repulsive,  and  the  outer  rings  are  again  attractive. 
Since  four  photon  bits  are  used  in  our  implementation,  and  since  each  ring  is 
either  attractive  or  repulsive,  two  rings  are  affected  by  a  simultaneous  parti- 


28 


tioning  of  the  space.  For  the  attractive  interaction,  there  are  5  type  of  orbits: 
bounce-back  with  A p  =  2,  clockwise  and  counter-clockwise  with  A p  =  1,  and 
clockwise  and  counter-clockwise  with  A p  =  \/ 3.  Consequently,  there  are  5  types 
of  repulsive  interactions,  which  are  just  the  conjugates  of  the  5  attractive  ones. 
Since  there  are  3  partition  directions  for  a  triangular  lattice,  it  takes  5*3=15 
partitions  of  the  space  to  completely  update  an  attractive  ring  and  a  repulsive 
ring  simultaneously.  To  compute  8  rings  therefore  takes  4*15=60  scans  of  the 
space.  Therefore,  since  the  local  collisions  require  a  single  scan,  it  takes  a  total 
of  61  scans  to  complete  one  time  step. 

A  two  dimensional  example  using  6  interaction  ranges,  with  an  underlying 
512  x  512  lattice,  of  this  time-dependent  crystallization  process  is  given  in  figure  9 
and  illustrates  the  type  of  molecular  dynamics  simulation  that  can  be  achieved 
with  a  more  complex  long-range  lattice-gas  algorithm.  The  resulting  crystal 
is  in  a  hexagonal-close-pack  configuration  since  we  have  strived  to  make  the 
coarse-grained  interatomic  potential  be  radially  symmetric.  This  long-range 
lattice-gas  model  had  six  interaction  ranges:  r  =  —2  ,— 7,  19,  21,  —24,  —26. 
Here  the  negative  sign  preceding  the  range  denotes  an  attractive  interaction  at 
that  range. 

12  Conclusion 

Although  the  lattice-gas  molecular  dynamics  algorithm,  described  above,  re¬ 
quires  61  scan,  which  is  quite  a  lot  of  scans,  this  implementation  only  requires 
10  bits  of  local  site  data.  Six  bits  are  used  for  the  momentum  states  and  4  bits 
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t=50 


t  =  500 


t  =  5000 


Figure  9:  Time  evolution  of  crystallization  in  a  2-dimensional  lattice-gas  with  multiple  fixed- 
range  2-body  interactions.  The  resulting  crystal  is  in  a  hexagonal-close-pack  configuration 
since  the  coarse-grained  interatomic  potential  is  radially  symmetric.  The  underlying  lattice  is 
512  x  512.  Started  with  a  uniformly  random  configuration  at  d  =  0.1.  Twelve  directions  are 
used  for  long-range  momentum  exchanges.  Grain  boundaries  and  defects  are  observed  during 
the  early  stages  of  the  crystal  formation. 
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for  messenger  states,  which  means  that  only  IK  byte  is  needed  to  store  a  long- 
range  rule.  Since  the  CAM-8  uses  a  16-bit  word,  there  still  remains  6  bits  of 
unused  local  data.  We  use  these  remaining  bits  to  hold  a  table  look-up  address. 
That  is,  since  the  size  of  the  CAM-8  look-up  table  static  random  access  memory 
(SRAM)  is  64K  bytes,  and  our  long-range  rule  only  requires  IK  bytes,  we  can 
store  up  to  64  long-range  rules  into  CAM-8’s  SRAM  memory.  Since  our  molecu¬ 
lar  dynamics  algorithm  decomposes  into  61  applications  of  the  long-range  rules, 
it  is  now  clear  why  we  have  chosen  to  use  up  to  eight  ranges.  Although  the  de¬ 
scription  of  our  implementation  may  sound  complicated,  in  fact  from  a  software 
development  point-of-view  it  was  the  most  direct  and  most  simple.  We  have 
traded  off  time  to  save  memory.  Yet  the  well-known  principle  of  computer  sci¬ 
ence  that  one  can  save  much  time  at  the  expense  of  using  more  memory  applies 
to  our  algorithm.  So  optimizations  of  our  algorithm  can  be  made,  particularly 
concerning  trading  off  an  increase  of  local  site  data  for  a  decrease  in  the  number 
of  needed  scans.  Clearly,  the  molecular  dynamics  algorithms  would  be  signifi¬ 
cantly  sped-up  if  it  were  implemented  say  on  a  64-bit  architecture.  Of  coarse 
in  this  case,  computation  by  table  look-up  would  be  inappropriate.  However 
by  making  use  of  lattice  isometries  and  rule  conjugation,  the  necessary  logic 
is  actually  quite  small,  as  evidenced  for  example  by  (37)  or  (52).  Therefore, 
a  programmable  logic  method  of  computation  should  be  better  than  the  table 
look-up  method  of  computation  currently  in  use  in  the  CAM-8  for  this  kind  of 
lattice-gas  algorithm. 
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